/*
 * Copyright 2003 Program of Computer Graphics, Cornell University
 *     580 Rhodes Hall
 *     Cornell University
 *     Ithaca NY 14853
 * Web: http://www.graphics.cornell.edu/
 * 
 * Not for commercial use. Do not redistribute without permission.
 */
#include "stdafx.h"

#include "Parallelogram.h"

PrimitiveList *Parallelogram::compile(int num,const GeometryObject **array,PoolAllocator *pool) const
{
//_asm int 3; //force breakpoint
  PackedParallelogramList *list = PackedParallelogramList::poolAlloc(num,pool);
  for(int k=0;k<num;k++) {
    assert(array[k]->getType()==GEOM_PARALLELOGRAM);
    const Parallelogram *t = (const Parallelogram*)array[k];
    list->setParallelogram(k,t->v0x,t->v0y,t->v0z,t->v1x,t->v1y,t->v1z,
      t->v2x,t->v2y,t->v2z,t->v3x,t->v3y,t->v3z,t->id());
  }
  return list;
}

void Parallelogram::expandBoxToInclude(AxisAlignedBox *box) const 
{
  box->addPoint(getV0());
  box->addPoint(getV1());
  box->addPoint(getV2());
  box->addPoint(getV3());
}

AxisAlignedBox Parallelogram::getBoundingBox() const 
{
  AxisAlignedBox retval(getV0());
  retval.addPoint(getV1());
  retval.addPoint(getV2());
  retval.addPoint(getV3());
  return retval;
}

bool Parallelogram::intersectsRay(const Ray &ray, IntersectionRecord *hitRec) const
{
  //first form edges
  F32vec4 e1 = getV1() - getV0();
  F32vec4 e2 = getV3() - getV0();
  //C is e1 X e2
  F32vec4 C = crossProduct(e1,e2);
  F32vec1 denom = rcp_nr(dotProduct(C,ray.direction));
  //S is vertex zero minus the ray origin
  F32vec4 S = getV0() - ray.origin;
  //t is the parametric parameter along the ray of the intersection
  //t = (S dot C) / (D dot C)
  F32vec1 t = dotProduct(S,C)*denom;
  //B is S X D
  F32vec4 B = crossProduct(S,ray.direction);
  //u is the first barycentric coordinate of the intersection
  //u = (e2 dot B) / (D dot C)
  F32vec1 u = dotProduct(e2,B)*denom;
  //v is the second barycentric coorinate of the intersection
  //v = -(e1 dot B) / (D dot C)
  F32vec1 negV = dotProduct(e1,B)*denom;
  //test to see if intersection is valid
  F32vec1 zero(0.0f);
  F32vec1 one(1.0f);
  F32vec1 negOne(-1.0f);
  F32vec1 test = cmpge(u,zero)&cmple(negV,zero)&cmple(u,one)&cmpge(negV,negOne)
    &cmpge(t,F32vec1(ray.tStart))&cmple(t,F32vec1(ray.tEnd));
  //I haven't thought of a better way of extracting mask from F32vec1 so we convert it to
  //a F32vec4, move the mask to an integer and then ignore everything except the low bit
  if (move_mask(F32vec4(test))&0x01) {
    Store_F32vec1(&hitRec->t,t);
    hitRec->primitiveId = id();
    Store_F32vec1(&hitRec->u,u);
    Store_F32vec1(&hitRec->v,zero-negV);
    hitRec->w = 0;
    return true;
  }
  return false;
}

bool Parallelogram::intersectsBox(const AxisAlignedBox &box) const
{
  //trivial accept if box contains any of the vertices
  if (box.containsPoint(getV3())
    ||box.containsPoint(getV2())
    ||box.containsPoint(getV1())
    ||box.containsPoint(getV0())) return true;
  //trivial reject if box does not intersect bounding box of triangle
  if(!getBoundingBox().intersects(box)) return false;
  //if simplistic tests were not definitive, we use the more expensive exact test
  Ray ray;
  IntersectionRecord hitRec;
  //get triangle edge vectors and take their cross product
  F32vec4 e1 = getV1() - getV0();
  F32vec4 e2 = getV3() - getV0();
  F32vec4 C = crossProduct(e1,e2);
  //now form the diagnonal of the bounding box that is closest to the plane normal
  F32vec4 posMask = cmpge(C,F32vec4Zero());
  F32vec4 S = Select(box.min,box.max,posMask);
  F32vec4 E = Select(box.max,box.min,posMask);
  //and test to see if this diagonal intersects the parallelogram
  ray.setAsSegment(S,E);
  if (intersectsRay(ray,&hitRec)) return true;
  //next check to see if any of the parallelogram's edges intersect the box
  ray.setAsSegment(getV3(),getV2());
  if (box.clipRay(&ray)) return true;
  ray.setAsSegment(getV2(),getV1());
  if (box.clipRay(&ray)) return true;
  ray.setAsSegment(getV1(),getV0());
  if (box.clipRay(&ray)) return true;
  ray.setAsSegment(getV3(),getV0());
  if (box.clipRay(&ray)) return true;
  //otherwise the box does not intersect this parallelogram
  return false;
}

//used to cache information about a parallelogram for fast repeated intersection in the next routine
class ParallelogramIntersectCache {
public:
  F32vec4 v0; //vertex zero
  F32vec4 e1; //edge 1
  F32vec4 e2; //edge 2
  F32vec4 C;  //e1 X e2
  F32vec4 invDC; //1.0 / (C*(max-min))   components used as denominator for axis-aligned edge vs. triangle intersection
};
//fast ray intersection routine for when some quantities can be precomputed and reused for 
//multiple queries against the same arallelogram.  We check a box edge encoded as a ray (actually line)
//for intersection with the triange.  It it does intersect, then we add the intersection point to
//the AxisAlignedBounding box.
static void boxEdgeIntersectsParallelogramFast(const Ray &ray, const ParallelogramIntersectCache &cached, int axis, AxisAlignedBox &box)
{
//  F32vec1 denom = rcp_nr(dotProduct(cached.C,ray.direction));
  F32vec1 denom = F32vec1(cached.invDC[axis]);
//  F32vec1 denom = F32vec1(1.0f)/dotProduct(C,ray.direction);
  //S is vertex zero minus the ray origin
  F32vec4 S = cached.v0 - ray.origin;
  //B is S X D
  F32vec4 B = crossProduct(S,ray.direction);
  //u is the first barycentric coordinate of the intersection
  //u = (e2 dot B) / (D dot C)
  F32vec1 u = dotProduct(cached.e2,B)*denom;
  //v is the second barycentric coorinate of the intersection
  //v = -(e1 dot B) / (D dot C)
  F32vec1 negV = dotProduct(cached.e1,B)*denom;
  //test to see if intersection is valid
  F32vec1 zero = F32vec1Zero();
  F32vec1 one(1.0f);
  F32vec1 negOne(-1.0f);
//  F32vec1 test = cmpge(u,zero)&cmple(negV,zero)&cmple(u-negV,one);
//  if (move_mask(F32vec4(test))&0x01) {
//  if (_mm_ucomilt_ss(test,zero)) {
  //I apologize for this convoluted expression.  We want to know if the result
  //of the compares is masks of ones (encoded as a -QNAN).  We could do this
  //quickly with an unordered compare, but Intel has not provided an intrinsic for
  //that instruction, so instead I have inverted the tests and used a greater-than-equal 
  //compare which will suceed when both of its arguments are zero and 
  //fail if either argument is a QNAN (such as a mask of all ones)
  if (_mm_ucomige_ss(cmplt(u,zero)|cmpgt(u,one),cmpgt(negV,zero)|cmplt(negV,negOne))) {
    //t is the parametric parameter along the ray of the intersection
    //t = (S dot C) / (D dot C)
    F32vec1 t = dotProduct(S,cached.C)*denom;
    F32vec4 pt = ray.origin + ray.direction*BroadcastFloat0(F32vec4(t));
    box.addPoint(pt);
  }
}

//test an edge of an AxisAlignedBox for intersection with parallelogram and if found,
//add intersection point to bounding box "bounds"
static void
testBoxEdgeP(AxisAlignedBox &bounds, const F32vec4 &start, const F32vec4 &end,
            const ParallelogramIntersectCache &cached, int axis)
{
  assert( ((start[0]==end[0])&&(0!=axis)) || ((start[0]<end[0])&&(0==axis)) );
  assert( ((start[1]==end[1])&&(1!=axis)) || ((start[1]<end[1])&&(1==axis)) );
  assert( ((start[2]==end[2])&&(2!=axis)) || ((start[2]<end[2])&&(2==axis)) );
  Ray ray;
  ray.setAsSegment(start,end);
  //we have already proved that the end points lie on opposite side of the parallelogram's plane
  //therefore the t value should always be between 0 and 1.  But for some reason, probably
  //roundoff error, it ends up slightly outside this range and is false rejected as not
  //hitting the triangles plane.  We avoid this by forcing an extended t range for the
  //test ray
  boxEdgeIntersectsParallelogramFast(ray,cached,axis,bounds);
}
//test all 12 edges of axis aligned box ("box") against the specified parallelogram 
//and adds any intersection points found to the bounding box ("bounds")
//This routine looks a little more complicated than it is because we have rewritten it
//to make extensive use of SSE simd instructions, but the idea is still the same
static void testBoxEdgesAgainstParallelogram(const Parallelogram *tri, const AxisAlignedBox &box, AxisAlignedBox &bounds)
{
  //precompute some values that we will use in intersection testing against the triangle
  ParallelogramIntersectCache cached;
//  {
    const Is32vec4 mask3Int(0,-1,-1,-1);
    const F32vec4 mask3 = _mm_castsi128_ps(mask3Int);
    F32vec4 v0 = tri->getV0();
    F32vec4 e1 = tri->getV1() - v0;
    F32vec4 e2 = tri->getV3() - v0;
    F32vec4 C = crossProduct(e1,e2);
    cached.v0 = v0&mask3;
    cached.e1 = e1&mask3;
    cached.e2 = e2&mask3;
    cached.C = C;
    //edges of box will be potential ray directions for intersection query
    //in each case only one component will be non-zero (they're axis aligned)
    //so we can precompute all values of denom = 1.0/dot(C,ray->direction)
    //in the components of 1.0/(C*(max-min)) and we don't have to compute
    //denom for each box-edge vs. triangle intersection query
    cached.invDC = rcp_nr(C*(box.max-box.min));
//  }
  const F32vec4 xMask = _mm_castsi128_ps(brite::VEC_0_0_0_N1);
  const F32vec4 yMask = _mm_castsi128_ps(brite::VEC_0_0_N1_0);
  const F32vec4 zMask = _mm_castsi128_ps(brite::VEC_0_N1_0_0);
  //use crossproduct to get plane equation and then test each of the boxes eight vertices
  //to see which side of the parallelogram's plane they lie on.  Then we only need to test
  //box edges whose endpoints lie on opposite sides of the plane
  F32vec4 minC = box.min*cached.C;
  F32vec4 maxC = box.max*cached.C;
  F32vec4 D = BroadcastFloat0(F32vec4(dotProduct(cached.C,cached.v0)));
  F32vec4 minmaxXminmaxY = unpack_low(minC,maxC);
  F32vec4 minmaxminmaxX = _mm_shuffle_ps(minmaxXminmaxY,minmaxXminmaxY,_MM_SHUFFLE(1,0,1,0));
  F32vec4 minminmaxmaxY = _mm_shuffle_ps(minmaxXminmaxY,minmaxXminmaxY,_MM_SHUFFLE(3,3,2,2));
  F32vec4 side0123 = BroadcastFloat2(minC) - D + minmaxminmaxX + minminmaxmaxY;
  F32vec4 side4567 = BroadcastFloat2(maxC) - D + minmaxminmaxX + minminmaxmaxY;
  //for each box edge that crosses the plane, test it for intersection with the parallelogram
  F32vec4 side0246 = _mm_shuffle_ps(side0123,side4567,_MM_SHUFFLE(2,0,2,0));
  F32vec4 side1357 = _mm_shuffle_ps(side0123,side4567,_MM_SHUFFLE(3,1,3,1));
  F32vec4 zero = F32vec4Zero();
  int testX = move_mask(cmple(side0246*side1357,zero)&cmpneq(side0246,side1357));
  F32vec4 side0145 = _mm_shuffle_ps(side0123,side4567,_MM_SHUFFLE(1,0,1,0));
  F32vec4 side2367 = _mm_shuffle_ps(side0123,side4567,_MM_SHUFFLE(3,2,3,2));
  int testY = move_mask(cmple(side0145*side2367,zero)&cmpneq(side0145,side2367));
  int testZ = move_mask(cmple(side0123*side4567,zero)&cmpneq(side0123,side4567));

  if (testX&0x01) {
    testBoxEdgeP(bounds,box.min,
                          Select(box.max,box.min,xMask),cached,0);
  }
  if (testX&0x02) {
    testBoxEdgeP(bounds,Select(box.max,box.min,yMask),
                           Select(box.min,box.max,zMask),cached,0);
  }
  if (testX&0x04) {
    testBoxEdgeP(bounds,Select(box.max,box.min,zMask),
                           Select(box.min,box.max,yMask),cached,0);
  }
  if (testX&0x08) {
    testBoxEdgeP(bounds,Select(box.min,box.max,xMask),
                           box.max,cached,0);
  }

  if (testY&0x01) {
    testBoxEdgeP(bounds,box.min,
                           Select(box.max,box.min,yMask),cached,1);
  }
  if (testY&0x02) {
    testBoxEdgeP(bounds,Select(box.max,box.min,xMask),
                           Select(box.min,box.max,zMask),cached,1);
  }
  if (testY&0x04) {
    testBoxEdgeP(bounds,Select(box.max,box.min,zMask),
                           Select(box.min,box.max,xMask),cached,1);
  }
  if (testY &0x08) {
    testBoxEdgeP(bounds,Select(box.min,box.max,yMask),
                           box.max,cached,1);
  }

  if (testZ&0x01) {
    testBoxEdgeP(bounds,box.min,
                           Select(box.max,box.min,zMask),cached,2);
  }
  if (testZ&0x02) {
    testBoxEdgeP(bounds,Select(box.max,box.min,xMask),
                           Select(box.min,box.max,yMask),cached,2);
  }
  if (testZ&0x04) {
    testBoxEdgeP(bounds,Select(box.max,box.min,yMask),
                           Select(box.min,box.max,xMask),cached,2);
  }
  if (testZ&0x08) {
    testBoxEdgeP(bounds,Select(box.min,box.max,zMask),
                           box.max,cached,2);
  }
}

AxisAlignedBox Parallelogram::estimateIntersectionBounds(const AxisAlignedBox &box) const
{
  AxisAlignedBox bounds = getBoundingBox();
  if (!box.contains(bounds)) {
    //if parallelogram is not contained in the specified bounding box, then we use a
    //more expensive method to compute the bounding box of the portion of the 
    //triangle that does lie within the specified bounding box
    Ray ray;
    bounds.resetToEmpty(); 
    //form bounding box of portion of triangle inside splits bounding box
    //by finding all possible extremal points and adding them to "bounds"
    //first we test for intersections of the triangles edges with the box
    ray.setAsSegment(getV3(),getV2());
    if (box.clipRay(&ray)) {
      bounds.addPoint(ray.getStartPoint());
      bounds.addPoint(ray.getEndPoint());
    }
    ray.setAsSegment(getV2(),getV1());
    if (box.clipRay(&ray)) {
      bounds.addPoint(ray.getStartPoint());
      bounds.addPoint(ray.getEndPoint());
    }
    ray.setAsSegment(getV3(),getV0());
    if (box.clipRay(&ray)) {
      bounds.addPoint(ray.getStartPoint());
      bounds.addPoint(ray.getEndPoint());
    }
    ray.setAsSegment(getV1(),getV0());
    if (box.clipRay(&ray)) {
      bounds.addPoint(ray.getStartPoint());
      bounds.addPoint(ray.getEndPoint());
    }
    //next test edges of the box for intersection with the parallelogram
    testBoxEdgesAgainstParallelogram(this,box,bounds);
    if (bounds.isEmpty()) {
      estimateIntersectionBounds(box);
    }
    assert(!bounds.isEmpty());
  }
  return bounds;
}

void PackedParallelogramList::setParallelogram(int index, float v0x, float v0y, float v0z, 
                                        float v1x, float v1y, float v1z,
                                        float v2x, float v2y, float v2z,
                                        float v3x, float v3y, float v3z, int id)
{
  assert((length > index)&&(index>=0));
  int i = index;
  paras[i].v0 = F32vec4(0,v0z,v0y,v0x);
  paras[i].e1 = F32vec4(0,v1z,v1y,v1x) - paras[i].v0;
  paras[i].e2 = F32vec4(0,v3z,v3y,v3x) - paras[i].v0;
  paras[i].id = id;
#ifdef PRECOMPUTE_PARALLELOGRAM_CROSSPRODUCT
  //C is e1 X e2
  paras[i].C = crossProduct(paras[i].e1,paras[i].e2);
#endif
}

void IntersectPackedParallelograms(class PrimitiveList *list, class Ray *ray, class IntersectionRecord *hitRec)
{
  F32vec1 tStart = F32vec1(ray->tStart);
  F32vec1 tEnd = F32vec1(ray->tEnd);
  {
    PackedParallelogram *paras = ((PackedParallelogramList*)list)->paras;
    for(int count = list->length;count>0;count-=1) {
#ifdef PRECOMPUTE_PARALLELOGRAM_CROSSPRODUCT
      F32vec4 C = paras->C;
#else
      //C is e1 X e2
      F32vec4 C = crossProduct(paras->e1,paras->e2);
#endif
//      F32vec4 denom = F32vec4(1.0f)/(ray->v_dirX*Cx + ray->v_dirY*Cy + ray->v_dirZ*Cz);
      F32vec1 denom = rcp_nr(dotProduct(ray->direction,C));
      //S is vertex zero minus the ray origin
      F32vec4 S = paras->v0 - ray->origin;
      //t is the parametric parameter along the ray of the intersection
      //t = (S dot C) / (D dot C)
      F32vec1 t = dotProduct(S,C)*denom;
      //B is S X D
      F32vec4 B = crossProduct(S,ray->direction);
      //u is the first barycentric coordinate of the intersection
      //u = (e2 dot B) / (D dot C)
      F32vec1 u = dotProduct(paras->e2,B)*denom;
      //v is the second barycentric coorinate of the intersection
      //v = -(e1 dot B) / (D dot C)
      //we compute negV = -v instead to avoid the extra operation
      F32vec1 negV = dotProduct(paras->e1,B)*denom;
      //test to see if any intersections are valid
      F32vec1 zero = F32vec1(0.0f);
      F32vec1 one = F32vec1(1.0f);
      F32vec1 negOne = F32vec1(-1.0f);
      F32vec1 valid = cmpge(u,zero) & cmple(negV,zero) & cmple(u,one) & cmpge(negV,negOne)
        & cmpgt(t,tStart) & cmplt(t,tEnd);
      int test = move_mask(F32vec4(valid))&0x01;
      if (test) { //there is a valid intersection
        if (hitRec->func) { //if we have a custom intersection responder, then use it
          (*hitRec->func)(ray,hitRec,F32vec4(t)[0],F32vec4(u)[0],-F32vec4(negV)[0],0,paras->id);
          tEnd = F32vec1(ray->tEnd); //reload in case it has changed
        } else { //otherwise the default is find first intersection
//_asm int 3; //force debugger breakpoint
          hitRec->t = ray->tEnd = F32vec4(t)[0];
          hitRec->u = F32vec4(u)[0];
          hitRec->v = -F32vec4(negV)[0];
          hitRec->w = 0.0f;
          hitRec->primitiveId = paras->id;
          tEnd = t;
        }
      }
      paras++;
    }
  }
}

